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We investigate the possibility of bistable lasing in microcavity lasers as opposed to bulk lasers. To that 
end, the dynamic behavior of a microlaser featuring two coupled, interacting modes is analytically investigated 
within the framework of a semiclassical laser model, suitable for a wide range of cavity shapes and mode 
geometries. Closed-form coupled mode equations are obtained for all classes of laser dynamics. We show that 
bistable operation is possible in all of these classes. In the simplest case (class-A lasers) bistability is shown 
to result from an interplay between coherent (population-pulsation) and incoherent (hole-burning) processes of 
mode interaction. We find that microcavities offer better conditions to achieve bistable lasing than bulk cavities, 
especially if the modes are not degenerate in frequency. This results from better matching of the spatial intensity 
distribution of microcavity modes. In more complicated laser models (class-B and class-C) bistability is shown 
to persist for modes even further apart in frequency than in the class-A case. The predictions of the coupled 
mode theory have been confirmed using numerical finite-difference time-domain calculations. 

PACS numbers: 42.65.Pc, 42.65.Sf, 42.55.Sa 
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I. INTRODUCTION 

In recent years, microlasers have been an object of growing 
interest in the photonics community because of a remarkable 
promise in both basic and applied research. Modern technol- 
ogy has facilitated fabrication of high-Q micro- and nanosized 
cavities (microresonators) in a vast variety of designs (mi- 
crodisks, -rings, -gears, -toroids, nanowires, nanoposts, and 
so on 1 1]). Lasers can be based on many of these set-ups as 
well as on different materials, e.g., semiconductors, impurity 
ions, or dye molecules. In addition, periodic nanostructures 
(photonic crystals, PhCs) can provide both cavity-based and 
distributed feedback resonators suitable for laser design M^i- 
The cavity size, which becomes so small as to be comparable 
to the operating wavelength, is what makes a microlaser phys- 
ically distinct from conventional ("bulk") cavities whose size 
is far larger. The small size limits the number of cavity modes 
that could take part in lasing, and at the same time greatly in- 
creases the influence of the cavity shape on the character of the 
modes. As a result, the mode structure becomes more com- 
plicated and heavily dependent on the specific cavity design. 
One is no longer able to describe the modes universally in 
an analytical manner. The variety of laser dynamics becomes 
much richer, which complicates the studies of microlasers to a 
considerable extent but at the same time can harbor interesting 
new effects. For example, one could look for new possibili- 
ties of bistable or multistable lasing |0], which would prove 
useful in many applications such as multiple- wavelength light 
sources, optical flip-flop devices or optical memory cells |5]. 

In the simplest case when two modes coexist in the same 
laser cavity (competing for the same saturable gain medium), 
three lasing regimes are usually considered t^. First, when 
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one of the modes has an advantage (e.g., larger Q -factor or 
better coupling to the gain), it simply dominates, becoming 
the only lasing mode (single mode lasing). Second, when the 
modes are well balanced (i.e., similar Q -factors and equally 
well coupled to the gain), they can both lase simultaneously. 
Such a coexistence can become possible because the modes 
with different frequency and/or spatial field pattern preferably 
interact with different gain centers. As a consequence, the 
spectral and spatial hole burning causes each mode to get sat- 
urated independently and allow the mode that happens to be 
weaker to catch up with the stronger one. Each mode saturates 
itself more readily as it does the other mode; in this sense, the 
coexisting modes are said to be weakly coupled {simultane- 
ous multimode lasing). Third, if the reverse is true, i.e., if 
each mode saturates the other mode before coming to its own 
saturation (the modes are strongly coupled), the weaker mode 
is quenched by the stronger one before it has any chance to 
catch up. Whichever mode has an initial advantage wins the 
competition and becomes the only lasing mode (bistable mul- 
timode lasing). The system can lase in either mode and is in 
this sense bistable. 

Trying to understand the physical origin of strong mode 
coupling brings about certain problems. It was pointed out 
from the beginning ^ that harmonic modes (such as longitu- 
dinal modes in bulk cavities) must always be coupled weakly 
because the antinodes of the field (the regions where the light- 
matter interaction is maximized) are spatially mismatched for 
different modes. Spatial hole burning would work similarly 
for any two modes with mismatched intensity distribution 
(such as transverse modes in bulk cavities). One of the ways 
to circumvent this limitation is to use degenerate modes with 
identical spatial intensity profiles, e.g., polarization degener- 
ate modes or counterpropagating modes in ring lasers. This 
can make the lasing bistable due to additional mode coupling 
through population pulsations iS, [SD- Alternatively, one can 
place a saturable absorber in addition to the saturable gain 
medium into the cavity |7, 9, 10]. Such an absorber can be 
naturally realized when only a part of the active medium is 
pumped. Both principles can be adapted for use in microlasers 



and are embodied in the form of polarization-bistable and ab- 
sorptive bistable laser diodes [11]. It has also been shown that 
two coupled lasers can achieve bistability if the output from 
each laser is directed to the other one and the feedback is re- 
duced to prevent formation of a compound cavity [12, JJ]. 
Later studies uA \1m give a detailed account on the stabil- 
ity and mode locking regimes of bulk coupled lasers based on 
nonlinear bifurcation analysis of the corresponding rate equa- 
tions. It is fundamentally problematic to achieve similar be- 
havior in microlasers where the modes share the same cav- 
ity. Recent achievements in the design of bistable multimode- 
interference laser diodes |16], though capable of bistable las- 
ing within a cavity of sub-millimeter size, still require sat- 
urable absorbers for the device to function properly. 

In the meantime, recent results show that there are yet unex- 
plored possibilities for bistable operation of microlasers. We 
have shown [17] that a cavity based on coupled defects in 
a PhC exhibits bistability without the need for saturable ab- 
sorption or similar additional mechanisms. The same idea 
was seen to work in lasers based on multimode nanopillar 
waveguides [18]. Similar results have been reported based on 
coupled microdisk |4] and coupled microring |5] resonators, 
the latter proposed for an ultrafast, ultralow-power optical 
memory cell design. Also, Ref. |19] reports that coupled 
multiple-feedback ring lasers can be brought to bistability by 
carefully selecting the feedback times, which may be more 
feasible in microlasers than the conventional gain-quenching 
scheme as in u3 . Finally, a time-independent multimode 
laser theory recently developed by H. Tiireci and co-workers 
11201 I21I1 reports that mode interaction can be very impor- 
tant in highly multimode nanostructure-based systems such 
as random lasers |22]. In view of this, there is a pronounced 
need to address the question of bistability in microresonators 
with their specific features such as complex cavity shapes 
and mathematically complex cavity modes taken into account 
consistently. Spatial hole burning should also be accounted 
for rigorously without reverting to averaging approximations, 
which are usually applied for coupled or semiconductor lasers 

In this paper, we consider the dynamics of two interact- 
ing modes in a microresonator-based laser. The semiclassical 
rate-equation model based on the Max well-B loch equations is 
used to model a laser-active medium. Coupled mode equa- 
tions are derived and analyzed for different classes of laser 
dynamics. Compared to existing accounts on mode dynam- 
ics and coupled lasers u^ [Tsl [23l U^, no specific form is 
assumed for either the cavity or the mode geometry. The 
spatial distribution of population inversion is taken into ac- 
count fully in terms of projections onto the modes' subspace 
(see |25]) for all classes of laser dynamics. The theory de- 
veloped here can be seen as complementary to the account in 
Refs. 1120, I21I] by being able to provide a description of time- 
dependent laser dynamics. Though they are rather different, 
both these approaches go beyond the third-order nonlinearity 
in the description of light-matter interaction. 

In the simplest case of class-A laser dynamics, the equa- 
tions suitable for analytical studies have been derived. As al- 
ready shown earlier for some particular cases (see, e.g., fl4\). 



we confirm that coherent mode interaction (population pul- 
sations) can result in bistable laser operation. We show that 
bistable lasing becomes increasingly more difficult to achieve 
as the intermode frequency spacing Acc; increases from zero. 
However, for microcavity modes with well-matched intensity- 
gain overlap the bistability window has been shown to be 
much greater (by up to several orders of magnitude with re- 
spect to Alj) than for harmonic bulk-cavity modes. A non- 
symmetric system, where one of the modes is given an ad- 
vantage through cavity design, is also investigated. We show 
that a parameter mismatch favoring one of the modes can be 
compensated for by an opposing mismatch in another param- 
eter that would favor the other mode. In the more compli- 
cated class-B or class-C cases, numerical studies of the ob- 
tained coupled mode equations have been carried out. The 
effects of increasing the pumping rate and/or Auj beyond the 
applicability limits of the class-A approximations are studied. 
Bistable lasing is seen to persist unless Aco becomes compara- 
ble to the width of the gain line. Even then, bistability can be 
further restored by increasing the pumping rate highly above 
threshold. The results obtained for the class-B/C microlaser 
systems in the framework of the coupled mode theory have 
been compared with full numerical finite difference time do- 
main (FDTD) calculations. At least for the system considered 
(coupled defects in a 2D photonic crystal as in Ref. iflTI] ). we 
demonstrate that the predictions of the theory are in a good 
agreement with the results of numerical simulations. 

The paper is structured as follows. In Section [III we derive 
the semiclassical coupled two-mode laser equations suitable 
for a wide range of microcavity modes. Only a few general 
assumptions about the cavity shape are made and no partic- 
ular form for the mode geometry is specified. The deriva- 
tion starts from the Max well-B loch equations and is carried 
out from the more general (class-C) through the intermedi- 
ate (class-B) to the most restrictive (class-A) laser dynamics. 
Specific issues pertaining to introducing the dynamics classes 
in multimode lasers are addressed along the way. The analy- 
sis of the equations obtained is then carried out in the reverse 
order. In SectionJIIl we analyze the class-A case, which, with 
some assumptions, turns out to be closely related to the stan- 
dard two-mode competition model |6]. The parameter win- 
dow of bistable operation is investigated in terms of the spa- 
tial and spectral mode properties. In Section [IVl class-B and 
class-C equations are numerically investigated, and the main 
differences with the class-A model as regards bistable lasing 
operation are discussed. Finally, Section |Vl summarizes the 
paper. 



II. COUPLED TWO-MODE LASER EQUATIONS 

A. Semiclassical laser equations and multimode expansion 

The semiclassical laser equations used in the present pa- 
per as a starting point are composed of three parts: (i) the 
laser rate equations, reduced to the equation for population 
inversion W of the laser transition; (ii) the equation of mo- 
tion for the macroscopic polarization density P of the laser 
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Figure 1 : (Color online) Schematic illustration of the mode frequen- 
cies (ooi,2 = oooT ^^) with respect to gain (Suj = ooa— (^o), as used 
throughout the paper. 



medium, obtained in a modified electronic oscillator model, 
and (iii) the scalar wave equation derived from the Maxwell 
equations. We consider two-dimensional (2D) systems, trans- 
lationally invariant in the z-direction, with TM light polar- 
ization, corresponding to a wide range of 2D photonic struc- 
tures. In this case the electric field is E(r) = Ez{x^ z)z, al- 
lowing us to restrict ourselves to the 2: -component of the field 
E{t^ t) = Ez{x^ y^ t). Applying the slowly varying envelope 
(SVE) approximation |6], the Max well-B loch system of equa- 
tions takes the form [i23|] 



^W{v,t) = 7||[i^-W(r,t)] + ^[^(r,t)P*(r,t)-^*(r,t)P(r,t)], 



(1) 



dt 



P(r,t) = -(7± + i(^)P(r,t)-^I^(r,t)^(r,t), 



(2) 
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Here W{r^ t) has the meaning of population inversion, which 
can vary spatially as opposed to Ref. [14] where it is assumed 
to be constant across the whole cavity. Further, R is the exter- 
nal pumping rate, /i is the dipole matrix element of the atomic 
laser transition, and the polarization and population inversion 
decay rates are given by 7^ and 7||, respectively. We con- 
sider a resonant system that features two eigenmodes with de- 
cay rates /t:i,2, phenomenologically accounted for by the pres- 
ence of a loss term n{v) in Eq. dS]). The mode frequencies are 
cc;i,2 = cjo T ^^. and the central frequency uoo is shifted with 
respect to the lasing transition frequency uOa by 5^^ = uOa — ^0 
with Acc;, 5^^ <C cjq, as shown in Fig. [T] We assume that the 
eigenmodes of the cold cavity have a spatial structure given 
by ^1,2(1*). The electric field E{y^ t) is then decomposed into 
the spatially dependent mode profiles 1^1,2(1*) multiplied by 
time dependent SVE functions £^1,2 (^) as 
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Here and further, ^± = ±iAujt. Following the approach in 
pS], we make a similar ansatz for the polarization, introduc- 
ing the amplitudes Pi,2(^) as 



P(r,t)e- 



-'^'' = [iii(r)Pi(t)e^+ +^i2(r)P2(t)e^-] e-'^'\ 

(5) 
The applicability of the expansion (O needs further justifi- 
cation. Eq. (|5]) assumes that polarization P(r, t) and the elec- 
tric field £^(r, t) have similar spatial profiles. This is strictly 
true only if the field intensity is small enough, e.g., if the 
pumping rate R is not very large. Otherwise, the polariza- 
tion gets influenced by the saturation terms that involve the 



population inversion VF(r, t), which itself cannot be spatially 
decomposed. These saturation terms would modify the spatial 
profile of P(r, t) outside the scope of Eq. (O. 

However, as Eq. (|5l) does not contain any explicit expan- 
sion in a series of nonlinearity orders with subsequent se- 
ries truncation, the constraint on the pumping rate R ap- 
pears to be much weaker than what is enforced by the usual 
near- threshold expansion [llfl |27l |28|], which explicitly re- 
tains only third-order nonlinearities in the hole burning in- 
teraction. In the extreme (single-mode) case, where Eq. ^ 
implies P(r, t) ex E{t, t) and thus carries the strongest ap- 
proximation, it can be shown that the coupled mode theory 
based on Eq. ([5]) leads to underestimation of the steady-state 
laser field intensity E{R). However, the character of the de- 
pendence E{R) is preserved for the values of R well outside 
the range of applicability of the near-threshold expansion (see 
||2Q|] ). Moreover, the dynamical behavior of the laser is also 
correctly predicted by the coupled mode theory employing the 
expansion J5l) both for one and for two modes (see our ear- 
lier work 1I25I] for a comparison with direct numerical simula- 
tions). 

That taken into account, in what follows we will use the ex- 
pansion (|5]), remembering that the results may deviate quanti- 
tatively and may be subjest to further checking as the pumping 
rate goes far above threshold. 



B. Class-C lasers 

In order to derive the equations for Ei{t) and Pi(t), one 
has to eliminate all the spatial dependencies from Eqs. ([T])- 



dS]). We begin by substituting Eqs. dU and (|5]) into Eq. dS]), 
assuming that the time dependence of the field envelopes are 
slow enough so that \dEj/dt\ <C ujj \Ej\. The modes Uj{t) 
are assumed to be orthonormal solutions of the homogeneous 
wave equation (c^V^ — e{r)u;'j)uj{r) = 0, which means that 
their integral across the cavity is 



c 



e{r)u;{r)uj{r) 



(6) 



As a result, the spatial derivatives in Eq. ^ can be eliminated. 
If e(r) = e is constant throughout the cavity (the bulk-cavity 
case 12311 ). the modes in Eq. ^ decouple rigorously, and one 
obtains 



dt ^ 



njEj 
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2eoe ' ' 



(7) 



This decoupling remains approximately true if the major part 
of the modes' energy are located in a material with the same 
dielectric constant, as is often the case in microcavities. For 
details we refer the reader to our earlier work [25]. A more 
complicated case of distributed feedback structures would re- 



quire additional spatial multiscale analysis, e.g., following the 
approach developed for photonic crystal lasers |24]. 

Eliminating spatial dependencies in Eq. ^ is simpler and 
requires substitution of Eqs. ©-(O with subsequent projec- 
tion onto the eigenmodes, i.e., integration j u*{. . .)d^Y over 
the gain medium: 



-Pi = - APi - i^ {EiWii + ^2^^126^^-) , 

dt a ^ ^ 






(8) 



P2P2 - i^ (^iT^2ie2*+ + E2W22) , 



where /?i,2 = (7± + i<5aj) ± ''^^'^ and Wy are the projections 
of the population inversion W(x,t) onto the corresponding 
modes 
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W,,(t) = e / d'vu'l{v)W{v,t)uj{v) (9) 

JG 

Analogously, by substituting Eqs. dU)-© into Eq. dB and ap- 
plying J ii*(. . .)iXjd^r, one can obtain the equations for Wij 
in the following form: 



J 



—Wij =-f\\{Rij 



W,. 



^ [E, {a]}P^ + a^jP^e'^-) + E, {a^P^^^' + o^^P^)] . 



(10) 



Here, Rij are related to R in the same way as Wij to VK(r, t), via Eq. dS- The coefficients a^^ are mode overlap integrals 



defined as: 






L 



d^YUl(Y)Uj(Y)ul^(Y)Un(v). 



(11) 



The integration in Eqs. dS and dHJ) is performed over the gain medium where e(r) = e is assumed to be constant. Apart from 
that assumption, the shape of the gain region itself can be arbitrary and does not have to be contiguous. The mode geometry can 
also be arbitrary unlike in the previous reports iHHI^^l, as the inter-mode and mode-gain overlaps are accounted for in terms 
of Qf^^ and Wij. Note that Eqs. dSl) and JTOl) with the definition ^ do not involve any approximations on the field or pump 
intensity beside the one associated with the validity of Eq. dSJ) as described above. Because of this, the full population inversion 
VF(r, t) cannot be written explicitly in terms of Wij{t) and 1^1,2(1*), in contrast to E and P, as in Eqs. dUl-dU). Also note that 
the rate equations dTOl) for the population inversion explicitly contain oscillatory terms, which originate from the beating in the 
superposition of the two modes with different frequencies uji and cc^2- 



C. Class-B lasers 



Equations dZj), dSl). and JTOl) govern the dynamics of the two spectrally close, interacting modes without any assumptions on the 
laser dynamics besides those needed for the SVE approximation. All these equations include a decay term with a characteristic 
decay rate for all the variables involved. The mode amplitudes Ej decay with the rate Kj associated with the Q-factors of the 
modes {Qj = ujj/ nj). The decay of all the population inversion projections Wij is governed by 7||. Finally, the polarization 
amplitudes Pj decay rates are complex, f3j = (7^ -\- iS^j) ^ lAco. This complexity directly results from the multimode character 
of the laser under study, and in the single-mode case Pj = 7^ . 



In the most general case of laser dynamics there are no restrictions on the decay rates Hij, 7|| . 7^ (so-called class-C lasers). In 
reality, however, the decay rates are governed by different physical processes and often belong to different time scales (class-B 
or class-A lasers, see 11150 ). which can make the analysis of the laser equations considerably simpler. 

Class-B lasers are defined by 7^ ^ 7|| ' ^i- ^^ ^^^ single-mode case, it would mean that the polarization relaxes and achieves 
saturation so fast that the polarization can be assumed to have no own dynamics and follows E and W adiabatically. 

In the two-mode case, where the polarization dynamics is influenced by the intermode spacing Acj, the introduction of 
the class-B approximations needs to be approached with greater care. Since the right-hand side of Eqs. ([8]) includes oscillatory 
terms on the time scale of 2Au;, one can eliminate the polarization only if these oscillations are much slower than the exponential 
decay due to 7^, i.e., 7^ ^ Acj. Note that this additional condition for class-B lasing, specific for multimode lasers, becomes 
especially important in microlasers where the small cavity size can place the modes much further apart from each other than in 
the bulk cavities. 

Under these assumptions, we can now eliminate the polarization adiabatically by assuming dPj/dt ^ 0. Hence, Eqs. ([8]) 
assume the form 
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which causes Eqs. d?]) to be modified as 
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Analogously, substituting Eq. (fT2l) into Eq. (flOl) one may obtain the equations for Wij . Since the population inversion W is real 
[see Eq. dB], it follows from Eq. ^ that W* 



dt 



^ that W*i 


= Wij , and in particular, W*j = Wjj . Hence, 
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Eqs. ([T3])-([T4l) are the governing equations for two-mode class-B lasers. Further knowledge about the modes in question can 
allow further simplification. A good example is the case when the modes are orthogonal not only withinin the whole cavity 
[Eq. ©I, but also in the gain region, e.g., if most of the cavity or at least the portion of the cavity with maximum mode energy 
is filled with the pumped gain medium: 
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<W^iW = Sij. 



(15) 



In this case the overlap integrals with one out-of-place index (a^^, ajl, etc.) will be negligible compared to the rest of the 
overlaps such as a^^, a**^, afj, or a^-. This allows to shorten Eq. ([T4b . which then assume different forms for symmetric Wjj 
vs. anti- symmetric projections Wij-Li'. 
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where i?i2 ^ Rjj due to the mode orthogonaUty and, as we remember, W21 = W12. Furthermore, if the modes are intensity- 
matched, i.e., assumed to have nearly equal intensity distribution in the gain region so that 
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then it follows from Eq. ^ that Wn = W22 = Ws and W12 = W^^ = Wa, as well as from Eq. ^ that a^ 



real, while a^- = a' can be complex. Hence, 
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D. Class-A lasers 

If one further assumes that (7^ ^)7|| ^ f^j (class- A lasers), the slowest- varying quantity becomes the mode decay. The 
population inversion follows the mode amplitudes Ej (t) instantaneously and can be eliminated, leaving us with only two equa- 
tions for the mode amplitudes. Similar to the way we have built the class-B approximation, the derivatives in Eqs. ([T4l) are 
dWij/dt ^ 0. In this case, Eqs. ([T6l)-([T7l) become 
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This is a system of linear algebraic equations that can be solved for Wij . We are aiming for equations with simple enough 
structure to be treated analytically, namely, equations for Ej with up to cubic-order non-linearity as analyzed, e.g., in ||^. 
Hence, we are looking for the solutions in the form 
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(23) 



neglecting terms with higher powers of E. Truncating higher-order nonlinearity corresponds physically to the case with low field 
intensities, i.e., just above the lasing threshold. Hence, at this point the near- threshold expansion is introduced as understood in 
numerous works 1 26l |27l [28] . We remark that this expansion is by far a stronger approximation than the one used in assuming 
the form ^ for the polarization. Hence, the class-B and class-C models described in the previous sections are valid for much 



stronger pumping, while the class-A description that follows is valid for pumping rates only slightly above threshold. Inserted 
into Eqs. (ITB-dJa. Eq. ^ yields 
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Note that the right-hand side of Eqs. (l2T1l-(l22b has terms of the form E^EnWij . Hence the same resuh could be obtained by 

/(fc) _ IT w(f=-i) with WJ--' 



^ ^ — up to W-^, ^ , as was done in 



solving the equation system Wij = L • Wij iteratively as W^j ^ = L 

Note that the presence of oscillatory exponents e^^^^ on the right-hand side of Eqs. ([T4l) . induced by beating of the field 
intensities, dictates that an adiabatic elimination can only be performed safely if 7|| :^ Auo. Unfortunately, this assumption is 
quite restrictive and makes the resulting class-A laser equations hardly applicable for any two-mode system beyond the case 
of spectrally overlapping modes unless the mode Q-factors become very high. However, Eq. (l24l) suggests that W12 should be 
oscillatory with frequency 2Auj. This is indeed the case, as confirmed by numerical solution of class-B or class-C equations. 
These oscillations (also called population pulsations) are the main reason why the condition dWi2/dt ?^ is valid only for 
vanishingly small Ac^;. By accounting for these pulsations explicitly, one can build class-A laser equations applicable for a wider 
range of Acc;. We introduce oscillatory terms e=^^^^^^ into 1^12: 



Wi2{t) = W2i{t) = Wa{t)e^'^+ + W:{t)e^'^- 



(25) 



where the envelope function Wa{t) supposedly varies more slowly than 2Acj and on the same time scale as Wjj (t). We can then 
reformulate the condition for adiabatic elimination of W12 in the form dWa/dt ^ 0. The algebraic equation for Wa analogous 
to (|22]) is then 
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(26) 

Note that unlike W12, Wa is explicitly complex due to the substitution 711-^ 7|| + 2iAw. Also note the disappearance of 
oscillatory exponents in Eq. (l26l) . compared to Eq. (l22b . Inserting Eq. (l25])-(l26b into (l22b and following the same near-threshold 
expansion as above, we obtain the final class-A equations 
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Eqs. ([27l) retain their 



applicability for a wide range of Acj up to Acj :^ 7|| and beyond. The only limitation is the requirement 7^ ;» Acc; needed to 
obtain the class-B equations. As was the case with the class-C to class-B transition, we see that the multimode case needs to be 
approached with care, since Auj represents an additional dynamical parameter (mode beating). It can play a significant part in 
laser dynamics and render some approximations invalid despite their validity in the single-mode case for the same parameters. 



III. BISTABILITY IN CLASS-A MICROLASERS 



A. Mode competition equations 



Now that the dynamics of a two-mode laser has been reduced to relatively simple class-A equations ([27]), the mode dynamics 
can be analyzed for possible steady-state and stable solutions. Eqs. dTTl) resemble the standard 2-mode competition equations 



(see li]): 
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Here, pj in the linear terms characterize the net unsaturated gain (minus cavity losses) for the mode j. The coefficients Ojj and 
Oij^i are self- and cross-saturation coefficients, respectively. These terms are fully similar in form and meaning to the widely 
studied case in [6]. The last terms, which are special to Eqs. (1281) , also contribute to cross-saturation but contain the phases of 
the modes, as well as an explicit oscillatory time dependence with frequency 4Acj. The expressions for all the coefficients can 
be obtained directly from Eqs. ([27]). 

Since Eqs. (l28l) include the phase of the modes explicitly, they can be separated into amplitude and phase equations. Substi- 
tuting Ej{t) = \Ej{t)\e''^^^^^ , one obtains: 



^ 1^1 1 = (Repi -Re^n |^i|' -Re^i2 |^2|') |^i| -Re {o[,e^'^^^-^^^e^^-) |^2|' |^i|, 
^ \E2\ = (Re p2 -Re ^21 |^i|' -Re^22 |^2|') |^2| - Re (^^le-^^^^^-^^^e^^^) \Eif |^2| , 



(29) 



dt 



^(^1 = (impi - Im^n |^i|' - lmO,2 |^2|') - Im (^^^e'^^^^-^^^e^^-) {E^f , 
^2 = (lmp2 - Im^2i l^il' - Im^22 |^2|') - Im (^0'^^e-^'^'^'-'^'^e^^+^ \Eif . 



(30) 



The amplitude equations (|29l ) now completely coincide in form with the usual two-mode competition [b] but contain the 
intermode phase difference A(p = (p2 — (fi sls 3. parameter and have the cross-saturation coefficients explicitly time-dependent. 
We can see that the amplitudes always achieve saturation due to a cubic non-linearity. The phase difference, however, may either 
become stationary, corresponding to phase-locked solutions, or be allowed to vary, in which case the solutions are said to be 
unlocked. 

In the limiting case of Acc; = one can show that there are two phase-locked solutions: one stable with Acp = 7r/2 and one 
unstable with Acp = 0. Without further assumptions as to the nature of the modes (such as those in some earlier works I!14l l26n). 
the general case is difficult to analyze due to explicit time dependence in the coefficients for non-zero Auj. In particular, Auj > 
causes Acp to undergo precession even in the locked regimes. As this precession becomes faster, one can no longer distinguish 
between locked and unlocked solutions. For sufficiently large Aco, the oscillations e^^^^^^ occur fast enough compared to the 
onset time scale, which primarily depends on k, rather than on Ac^;. In this case the modes appear always unlocked (mentioned 
in [26] as a "natural tendency" for different- frequency modes), and the effects of the phase terms can be averaged out. Our 
numerical estimations show that this is possible if Auj > 10~^a>:. The case Auj <C k., corresponding to spectrally overlapping 
modes, is outside the scope of the present paper anyway as there can be additional channels of mode coupling (e.g., the Petermann 
excess noise ll29ll ). Thus, we will henceforth ignore the phase terms in Eqs. (|28])-(|301) and rewrite Eq. (|27]) as 
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B. Conditions for bistable lasing: Mode coupling 



With the phase terms dropped, Eqs. (|3T1) represented in the 
amplitude form analogous to Eqs. (|29b can be analyzed fol- 



lowing the standard procedure fu]. The primary parameter 
that determines the nature of mode competition is the mode 
coupling constant 



C = Rei9i2Rei92i/Rei9iiRei922, 



(32) 
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Figure 2: (Color online) The dependence of Co; in Eq. (l35l l on 
7ll and Ao;. The dashed lines are the isolines for Cc^ = 1 and 
Cuo — 9/4. The dotted lines approximately mark the applicability 
limits of class- A equations. 



which is the ratio of cross-saturation and self- saturation coef- 
ficients. It is commonly known that the cases of simultane- 
ous two-mode lasing and bistable lasing are characterized by 
C <\ (weak mode coupling) and C > \ (strong mode cou- 
pling), respectively |6, 7]. Assuming that the pumping does 
not favor either of the modes so that Ri = R2 = R, as well 
as cji ?^ CJ2 = ^ ^ Acj, we can substitute the explicit form 
of the coefficients from Eq. ([3T]) into Eq. ([32l) . As a result, we 
have found that C can be factored as 



G — Cq;GcJ 



(33) 



The first factor C^, which originates in the spatial hole 
burning, has the form 



^a 



^2 



aiia22 



(34) 



In the simplest case when the modes are intensity matched as 
in Eq. ([TS]) so that all aij = a, it follows that Cq, = 1. Oth- 
erwise, it can be proven that C^ < 1. The second factor C^, 
which results from population pulsations and becomes identi- 
caly unity if those pulsations are neglected, has the form 



4Acj2 



a 






0,| 



(35) 



The dependence of Cu; is shown in Fig. [2l We can see 
that Co; < 4 for Au <C 7|| and C^j ^ I for 7|| < Acc; <C 7^. 



The transition between two limiting cases (Cuj = 1 and 
Coo = 4) occurs rapidly around Acc; ^ 7|| . Note that as Aco in- 
creases, Cou approaches unity from below, so there is a critical 
value Acj'^'^^ ^ a/7±7||/2 for which C^j = I. Hence, in the 
ideal case of intensity matched modes [Eq. ([T8l) 1 bistability is 
possible for Acc; all the way up to Au;^'^\ The limiting case 
of C = 4 is known to be realized for the ideal case of coun- 
terpropagating modes in ring lasers or modes with orthogo- 
nal polarizations, which are fully intensity matched and have 
Acj ?^ |6]. 

If, however, the modes are considerably mismatched, then 
Cuj must be significantly larger than one to compensate for 
a small Ca and thus keep the overall mode coupling con- 
stant above unity to achieve bistable lasing. For example, it 
can be shown that ID harmonic (e.g., longitudinal) modes al- 
ways have Ca = 4/9 for different frequencies. This means 
that the line of critical values for Aoo^^^^'^ ^ 7||/2 up to 
where bistability is possible lies much deeper than the line 
of Acj*^^^ (see Fig. O. Taking into account that the fre- 
quency shift between longitudinal modes is related to the cav- 
ity length as Acc;(buik) = ttc/L, one easily obtains the "rule of 
thumb" for minimum cavity length of a ID bistable bulk laser: 
Lmin — 27rc/7|| . For realistic laser media, Lmm is found to be 
prohibitively large, from around 2 m for semiconductors and 
up to 200-300 km for Nd:YAG |30]. This explains why it is 
so difficult to achieve bistable lasing for different-frequency 
modes in a bulk cavity: unless the cavity is extraordinary big, 
Auj is large enough to bring C^j so close to unity that any in- 
tensity mismatch causes Cq, < 1 and brings the laser back 
into the weak-coupling (simultaneous lasing) regime. The 
only notable exception is the case when the modes are quasi- 
degenerate with Acc; ^ 0, such as counterpropagating modes 
in ring lasers or modes with orthogonal polarization, and it 
is in these special cases that bistability could indeed be ob- 
served. 

In a microcavity, however, the modes can be made very 
nearly intensity matched by a carefully chosen resonator de- 
sign (e.g., coupled cavity-based, see |17]). In addition, many 
designs allow to control the frequency separation between the 
modes more or less independently from other model parame- 
ters. This opens up a whole new frequency range Acj*^^/^^ < 
Aco < Acc;*^^^ available for bistable laser design, which can 
encompass several orders of magnitude for Auj (see Fig. O. 
This range becomes available in microlasers because the pos- 
sibility to bring the modes to intensity matching is far greater 
than in bulk cavities, owing to a greater variety of cavity 
shapes and a more complicated nature of the modes involved. 

Finally, from Eqs. ([3T]) one can also see the physical mecha- 
nism of bistable lasing in the class-A case. It is due to the (os- 
cillatory) component W12 that there is an addition to the cross- 
saturation coefficients Oij^i. Without this addition, C would 
simply coincide with Ca and all possibility for bistable opera- 
tion would be excluded. Hence, it is the coherent mode inter- 
action effects such as population pulsations or four- wave mix- 
ing iTj] that make bistability possible. Incoherent effects (e.g., 
spatial hole burning, which is only manifest in Ca) can either 
allow or suppress it. As a result, an interplay between coher- 
ent and incoherent mode interaction processes is employed to 
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achieve bistable microlaser operation. 

As an example, we have plotted the dynamics of mode 
amplitudes Ej (t) as a numerical solution of Eqs. ([311 for 
bulk-cavity (Ca = 4/9) vs. coupled-cavity (Ca = 0.9) 
modes (Fig. O. Also shown are the temporal flow diagrams 
(i.e., projections of the solutions onto the |£^i| vs. \E2\ 
plane for different initial conditions of the cavity (the ratio 
e[^1 = \Ei{0)\ : |£^2(0)|). AH other parameters are kept 
constant, as described in the caption. If the modes are mis- 
matched (Fig.[3t) and C < 1, the laser saturates to the two- 
mode simultaneous lasing (|^i| = \E2\ = const) regardless 
of the initial conditions. Only this fixed point is stable. How- 
ever, if the modes are well matched (Fig. [3j)) so that C > 1, 
the laser saturates to a single-mode lasing as the initially 
stronger mode quenches its weaker counterpart and becomes 
dominant. There are two stable fixed points on the diagram: 
1^1 1 = const, \E2\ = and |£;i| = 0, \E2\ = const. The 
previously stable fixed point becomes unstable, and the line 
1^1 1 = 1^2! marks the separatrix between the stable points' 
domains of attraction. The mode that has an advantage in 
the beginning determines the domain of attraction for the sys- 
tem and hence the fixed point the system will converge to, as 
the separatrix cannot be transcended without an external influ- 
ence. These examples show that bistable lasing is possible in 
microlasers in such cases where only two-mode simultaneous 
lasing can be observed for bulk-cavity modes. 



C. Conditions for bistable lasing: Mode mismatch 

Up to now, we assumed that none of the modes is favored 
either by the cavity or by the gain, i.e., ki = K2, Ri = R2, 
ail = 0^12, and S^ = 0. In this case, as seen in Fig. [S] the 
two-mode lasing fixed point (labeled FP2), whether stable or 
unstable, is characterized by |£^i| = \E2\. This means, on 
the one hand, that in the simultaneous-lasing case both modes 
lase with equal intensity (Fig.[3t), and on the other hand, that 
in the bistable regime even a slight edge given to either mode 
in terms of initial conditions will bring this mode to lase. It 
is equally e asy to "select" or "switch" either mode by locking 
into it yj, 113]. This is illustrated in Fig. [3h by the fact that 
each stable fixed point has an equally large domain of attrac- 
tion. 

In a more general case, the ratio between mode intensities 
at FP2 Ii:2 = \Ei\ / I £^2 1 will change to reflect an advantage 
given to either of the modes. For example, even a slight mis- 
match in the mode Q-factors causes /i:2 to deviate from unity 
(Fig. [4]). Similar to the explanation given above, this may 
mean two things. In the simultaneous-lasing case, it simply 
means that once the laser achieves saturation, one mode has a 
greater amplitude than the other, e.g., |£^i| > |£^2|for/i:2 > 1 
(Fig.|4^). In the bistable case, it means that the domains of at- 
traction for the two modes change their size in phase space 
(Fig.llJ)). If for example /i:2 < 1, then the domain of attrac- 
tion for Mode 1 becomes larger, so Mode 1 is "in favor" as 
a result. In the bistable regime, a shifted FP2 means that the 
mode with a smaller domain of attraction is out of favor and 
thus harder to bring to lasing. For example, if FP2 is placed 



symmetrically, initial mode amplitude ratios E[.2 of 3:2 and 
2:3 bring the first and the second mode to lasing, respectively 
(Fig.[3]3). For asymmetrically placed FP2, the same two cases 
for initial condition both result in the lasing of the first mode 
(Fig.|4j)). To be able to target the smaller domain, one has to 
excite the out-of- favor mode exclusively, which might be dif- 
ficult experimentally. Hence we will further aim at finding the 
manifold of the system parameters for which I1.2 = 1. 

Whenever C 7^ 1, the general expression for the mode in- 
tensity ratio at FP2 /i:2 can be written as |6] 



/l:2 



Re piRe 6^22 - Re p2Re O12 
Re p2Re On — Re piRe O21 



(36) 



By substituting the coefficients in Eqs. ([3T1) one can obtain 
an explicit analytic expression for /i:2. Unfortunately, this 
general expression is very bulky and we will first investigate 
its behavior in several simplified cases. Let us introduce the 
perturbations in the form 



1^1,2 = /^(l ± ^/^), Rl,2 = ^(1 ± Sa), 



(37) 



from where it follows [see Eqs. ([9]) and ([TT1) 1 that a 11,22 

a(l ± Ja)^. Now if ^o; = ^a = 0, Jh: 7^ 0, /i:2 is givcn by: 



\-2 — 7 r • 



(38) 



Likewise if 5^^ = S^ = 0, Sa 7^ 0, then the expression is 
somewhat more complicated and reads: 



(a) _ (gg + Ca6l + CgS^) + {bg + dgSl) Sg 
^•^ {ag + CgSl + CgS^) " {bg + (IgSl) Sg ' 

Finally, if if Sg = 6,^ = 0, S^ 7^ 0, then 



(39) 



^(a;) ^ (gg; + c^Sl + e^S^) + {b^ + d^6l + US^) S^ 



(40) 

The coefficients in Eqs. (|38])-(|4Q1) are complicated polynomial 
functions of the dynamical parameters 7^ and 7||, the inter- 
mode frequency separation Acj, the measure of mode inten- 
sity mismatch v = ai2/a (which ranges from to a maxi- 
mum value of 1 — (5^ so that Cg < 1), and the pumping rate 
normalized to the threshold pumping R±y = 2Rgu;/{'yj_K,) 
1 32]. Note that k. itself does not enter these equations explic- 
itly. It does, however, impose a limitation Auj > 10~^/^ so 
that the phase terms in Eqs. (l27l) can be averaged out. 

From the structure of Eqs. ([38l)-(l4Qb one can see that /i:2 = 
1 for Sg = S^ = Su; = 0, as should be expected. If any one 
of the perturbation parameters (S^^^^g collectively referred to 
as 6) is non-zero, /i:2 deviates from unity. Obviously, chang- 
ing the sign of all non-zero S causes /i:2 -^ l/h:2- If favoring 
one of the modes (by any means) results in a certain asymme- 
try in lasing quantified through a non-unity /i:2, then favor- 
ing the other mode in the same way and by the same amount 
naturally causes the same asymmetry with respect to the other 
mode 1 33]. This suggests that one can choose more than one S 
to be non-zero in such a way that the shifts of FP2 caused by 
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Figure 3: (Color online) Flow diagrams of numerical solutions of Eqs. dSTT) in the \Ei\^ vs. |^2|^ plane for different initial mode ratios 
E[.2 — 1^1 (0)1 : 1^2(0)1 for (a) poorly matched modes {Coc — 4/9) and (b) well-matched modes {Ca — 0.9). The other parameters are the 
same (Ao; = 7|| = 10~^7±, ^1:1 = ^1:2 = 10~^7±), chosen such that 1 < Coo < 9/4. Solid circles and squares denote stable and unstable 
fixed points, respectively. The thin dashed line denotes the separatrix in the bistable lasing case. The right panels show the example mode 
dynamics |£^i,2 (t) | for two chosen values of E{.l slightly in favor of each mode. 



individual perturbations would cancel each other out. As a re- 
sult, one could achieve the resulting /i:2 equal to or close to 
unity, and the restrictions on the initial conditions would be 
lifted. 

Fig. [5] shows the manifold of the points /i:2 = 1 in the 
3D perturbation space {^uj^^k^^^o) for different parameters 
as a numerical solution of Eq. ([36l) . We can see that this 
manifold is an open surface. Hence, if a mismatch in one 
respect is unavoidable, it can be compensated for by engi- 
neering the other two perturbation parameters. Note that in 
the {5^] 5 ex) plane the mismatch compensation {1 1-2 — 1) is 
achieved when 5^ ?^ (5q, . This is easily understood if one re- 
members that the linear terms in Eqs. ([27]) have the structure 



ft ^ C^j 



C,R(\±^a) - n(\±^^). On the other hand. 



in the {^^\ (5^) plane, compensation is generally achieved for 
the oppositely-signed (5^ and 5^^. This is in agreement with an 



intuitive guess that, e.g., 5,^ < {k\ < ni) and (5^; < (the 
gain frequency cj^ < ^0 is closer to cji than to CJ2, see Fig.[T]) 
both give an edge to the first mode, so oppositely-signed b are 
needed to maintain balance. However, in the vicinity of the 
origin the surface can be folded, so that it crosses the origin 
with the opposite slope and compensation is achieved when ^^ 
and (5,^ have the same sign. Since perturbations ^o;, ^a, and (5,^ 
can have different physical origin and can be varied more or 
less independently by a proper choice of a gain medium and 
a cavity configuration, one can deliberately engineer a micro- 
laser to achieve bistable operation even if the idealized, unper- 
turbed case is difficult to realize experimentally. An example 
of such compensation is changing the mode frequencies with 
respect to gain (which can be done straightforwardly just by 
scaling the cavity) to help offset the difference in mode Q- 
factors, as shown numerically in our earlier work [18]. 
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Figure 4: (Color online) Same as Fig.[3]for mismatched mode decay rates (ki < K2, S^ 
arrows show the shift of the two-mode fixed point FP2. 



-0.42%) in favor of the first mode. The dotted 



To achieve a fully symmetric placement of FP2, one needs 
to bring three perturbation parameters into a relation. Be- 
cause all these parameters show only an indirect dependency 
on the cavity design and/or gain medium choice, the precise 
control of them may still be a challenging task. Hence, it is 
worthwhile to investigate to what extent the relations for ideal 
compensation can be violated so that bistable operation is still 
possible (albeit, as shown above, at the cost of stricter require- 
ments on the initial conditions). In terms of Fig. [51 that means 
how far one can deviate from the /i:2 = 1 surface and still 
lase into either of the modes on demand. 

From Eqs. (I38l)-(l40b one sees that a sufficiently high value 
of any 5 will cause either the numerator or the denominator 
in /i:2 to approach zero. On the flow diagram, this corre- 
sponds to the FP2 meeting the coordinate axes. Increasing S 
further causes /i:2 to become negative. The FP2 vanishes and 
the system finds itself in the single mode lasing regime (see 
||6t]). That sets an upper limit for any \S\ beyond which no 
bistable lasing is possible any more. 



More generally, the domain in space {Soo;S^;Sa) 
where /i:2 > comprises the possible perturbation parameter 
window where both modes can lase (either simultaneously or 
subject to bistability-induced switching, as depends on C). 
This domain, called the FP2 existence domain, is shown in 
Figs. [6]-[7l The existence domain, bounded by the surfaces 
defined by /i:2 = and /i:2 = oo is seen to surround the 
"perfect matching" surface /i:2 = 1. The domain boundaries 
appear to slide inwards as the pumping rate increases (Fig. [6]), 
which enlarges the FP2 existence domain around the point 
6 = 0. Also, the domain shrinks rapidly as the boundaries 
close around /i:2 = 1 when C approaches unity (Fig.|7]). The 
latter can be intuitively understood because C ^ 1 represents 
a delicately balanced system, so that even a slight mismatch 
is enough to throw the system heavily out of balance. Such a 
property is clearly a misfortune for the microcavity- specific 
bistability range reported above, as it relies on the situation 
when Coo exceeds unity only slightly. However, increasing 
the pumping appears to counteract this disadvantage, at least 
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Figure 5: (Color online) Manifolds of points /i:2 = 1 in the 3D 
perturbation space (^c^; ^k; fc) for (a) C = 1.10 > 1, (b) C = 
2.06 ~ 2, and (c) C = 3.84 < 4. The four surfaces correspond 
to four values of the pumping rate (R/Rthr = 1, 1-25, 1.5, 1.75), as 
indicated in the panel (b). 



Figure 6: (Color online) Boundaries of the FP2 existence domain 
1 1:2 > (dark gray) and the /i:2 = 1 surface lying inside that 
domain (light green) for C < 4 (as in Fig. |5j:): (a) at threshold 
(i?/i?thr = 1), (b) 10% above threshold (R/Rthr = 1.2), and (c) 
20% above threshold (R/Rthr = 1.2). 



for smaller S^j (see Fig. [6]). We believe that it is this effect 
that enabled us to observe bistability in earlier numerical 
simulations |17|, ilSj involving the laser operating highly 
above threshold. 

The practical conclusion to this section is that there are two 
theoretical requirements needed to achieve bistable lasing. In 
the first place, FP2 needs to exist on the flow diagram, as im- 
posed by /i:2 > 0. In the second place, once FP2 exists, the 
mode coupling constant must exceed unity (C > 1), as dis- 



cussed before. First (Sec. lIIIBI) . we have shown that in com- 
parison to bulk-cavity lasers microlasers exhibit a much wider 
parameter window characterized by C > 1, because the mi- 
crocavity modes can better fulfill the intensity matching con- 
dition ([TSl) . Secondly (Sec. lIIICl) , we have shown that there is 
an extended domain in the 3D perturbation space {S^;S^; 6a) 
where /i:2 > 0. Inside this domain, the closer /1.2 is to unity, 
the easier it is to realize bistability-based laser mode switch- 
ing experimentally. We have shown that /i:2 can be brought 
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Figure 7: (Color online) Same as Fig.[6]for constant R at 10% above 
threshold and (a) C > 1, (b) C :^ 2, (c) C < 4, corresponding to 
the cases in Fig.|5t-c. 



close to 1 by choosing a combination of perturbation param- 
eters that would compensate each other's advantage given to 
either mode. 



IV. CLASS-B/C MICROLASERS 

The elegance of the class- A case considered in the previous 
section is that Eqs. (IZTl) can be subject to analytical investiga- 
tion based on a comparison with Eqs. (l28l) |6]. Once a laser 
with a more complicated dynamics needs to be examined, 



more complicated systems of equations [six equations ([T3])- 
([T4l) for class-B or eight equations dTt-ifTOl) for class-C] need 
to be dealt with. Although attempts at analytical investigation 
of class-B equations are known (e.g., a near-threshold expan- 
sion of population inversion as proposed in |26]), only numer- 
ical solution seems to be applicable in the general case when 
no specific assumptions on the cavity or mode geometry are 
implied. Since all the equations are ordinary differential, such 
a numerical solution can be carried out with relative ease - the 
computational demands are far lower than a direct numerical 
integration of the Max well-B loch equations by means of an 
FDTD-likescheme|25]. 

A systematic investigation of class-B/C microlasers would 
be too lengthy to include in the present paper and will be the 
subject of a forthcoming publication. In this section we will 
outline the main differences in the behavior of such lasers 
compared to the previously studied class-A case as regards 
bistable lasing. 

We begin with a comparison of the laser classes in the near- 
threshold regime. As should be expected, the solutions for 
all classes display full coincidence if the class-A approxima- 
tion 7x ^ 7|| ^ ^ holds (note that this condition is rather 
restrictive in microlasers, requiring a careful choice of the 
gain medium as well as the cavity design). The mode dy- 
namics Ej{t) start to exhibit differences whenever 7|| or k 
are increased out of the class-A approximation. The differ- 
ences, however, are relatively minor, manifesting themselves 
mainly in the character of the transition process. In most 
cases, the mode coupling constant C as defined for the class-A 
in Eqs. ([33l)-([35l) continues to predict the laser dynamics cor- 
rectly (C < 1: simultaneous lasing, C > 1: bistability) even 
outside its strict range of applicability, although the behavior 
of Ej (t) can be quite different during the transition period. 

As discussed above, the class-B equations ([T3l)-([T4b do 
not involve a near-threshold approximation, it becomes pos- 
sible to consider a greater range of pumping rates, including 
regimes far above threshold, which are often left out of the 
picture in a construction of a multimode laser model ll23n . 
Comparison of the numerical results for class-B vs. class- 
C equations show that as long as the class-B prerequisites 
7± ^ 7|| , f^j hold, the results are similar, unless the condi- 
tion 7^ ;» Aco is violated. This agrees well with the ear- 
lier discussions in Section HTCl The differences appear not to 
be qualitative, but quantitative only, manifesting in the exact 
shape of the \Ej (t) \ dependence. The overall outcome of the 
mode interaction largely remains the same. To summarize, the 
main effect of the class-A to class-B transition in the context 
of studying bistable lasing is the inclusion of larger pumping 
rates R, while the main effect of the class-B to class-C transi- 
tion is the inclusion of larger frequency mode separations Aco. 

The increase of the pumping rate in a class-B laser is 
known to change the saturation character of the mode am- 
plitudes. The non-instantaneous relaxation of the population 
inversion with respect to the cavity field gives rise to spiking 
(for smaller R) or relaxation oscillations (for greater R) in the 
dependence Ej {t). A still stronger pumping (several orders of 
magnitude above threshold) causes the oscillations to vanish, 
as reported in an earlier work [ 25.1 . 
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Figure 8: (Color online) The dependence of lasing regime on pump- 
ing rate R and intermode frequency separation Aco in a class-C 
laser model. The parameters are Ai:j ^ 0.l7^ and 7|| ^ 10~^7^, 
as used for numerical simulation in |17]. The density plot shows 
the quantity_||^i_(t)| - |^2(t)|| /max (|^i(t)| , |^2(t)|) for large 
t::>7^^,77^,Ai:~^. Near-zero (light) values indicate two-mode (si- 
multaneous) lasing while near-unity (dark) values indicate one-mode 
(bistable) lasing. The lasing threshold depending on Ao; is marked 
with the dotted line. Numerical results of the FDTD simulations for 
coupled-defect structures (Fig. [9]) are superimposed over the density 
plot. Circles (red) and squares (yellow) show the location of points 
where simultaneous and bistable lasing, respectively, was observed 
in the mode dynamics during simulations. 



More interestingly, an increase of R can restore bistable 
lasing in the cases when simultaneous lasing is observed just 
above threshold. Fig. [2] suggests that there should be no bista- 
bility in the area around Aco c^ j± . The numerical solution 
of the class-C equations shows that this is indeed the case for 
smaller R. However, if the pumping is increased beyond a 
certain critical value R^ a transition from simultaneous to 
bistable lasing occurs (Fig. [8]). This effect was reported earlier 
1I25I1 with the observation that bistability ensues when pump- 
ing becomes so large that relaxation oscillations disappear. 
Our further investigations have revealed that this observation 
was rather a coincidence, and Re scales with Aco (Fig. [8]), 
bifurcating from threshold at approximately the point where 
Cuj = I according to Eq. (|35]) . This falls in line with the 
result of the previous section that a stronger pumping is capa- 
ble of restoring bistability where it has been deteriorated by 
adverse effects of insufficient mode matching. 

Because applicability of the expansion (|5]) and sometimes 
even of the SVEA f3V\ may become questionable far above 
threshold, we have carried out a comparison of Class-B/C 
results with direct numerical simulations. As previously de- 
scribed in Ref. [.25.1 . a space-time FDTD solver was coupled 
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Figure 9: (Color online) The family of structures used in numerical 
FDTD simulations, based on two coupled defects in a 2D photonic 
crystal lattice |17, 25]. Placing a different number of lattice rows 
between the defects (1-5), the intermode frequency separation Aco 
can be changed. 



to the four-level laser rate equations in order to model the re- 
sponse of a laser medium. A 2D photonic crystal lattice with 
two coupled defects um was used as a model system (Fig. [9]). 
Both defects are filled with four-level gain medium and con- 
tain a dipole source in the centre . By exciting these sources 
with varying amplitude/phase relations, the two modes (sym- 
metric and antisymmetric |25]) can be excited in any pro- 
portion and thus the initial state of the resonator can be con- 
trolled. By changing the number of lattice rows between the 
defects from 1 to 5, one can change Aco from ~ 7^ down to 
~ 10~^7x . The waveguides coupled to the defects form the 
primary channel for the radiation to leak out of the resonator. 
Care was taken that the mode Q -factors remain approximately 
the same across the whole family of structures. 

The results of the FDTD simulation runs are superimposed 
in the phase diagram in Fig.[8l For all values of Acc;, the tran- 
sition between simultaneous and bistable lasing was found ap- 
proximately around Re as predicted by the analytical theory. 
For larger Aco the correspondence is better because smaller 
Aco and R require much longer times to get to the steady state 
and there is an increased sensitivity to mode mismatch (see 
Fig. [6]). Hence it becomes more difficult to establish the transi- 
tion point between simultaneous and bistable lasing with good 
accuracy. 

In Figs. [To] and [TT] temporal laser dynamics in numerical 
simulations and the Class-C model are compared. We ana- 
lyze the electric field in the center of either defect Vc. For 
FDTD, it is monitored directly by recording the field at the 
corresponding point in space E{Yc-,t). To reduce the ex- 
cessive amount of data, we sample the field only at the lo- 
cal maxima, so that an envelope over the light oscillations is 
plotted. In the case of the coupled mode theory, the same 
quantity is obtained from the mode amplitudes ^1,2 (^) using 
Eq. dH) as Er{t) = u^^^ {Ei{t)e-''^'^ + E2{t)e-''^^^) where 
^max = max [ii(r)], assuming the modes are normalized ac- 
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Figure 10: (Color online) Comparison between laser field dynamics obtained by (a) the Class-C equations and (b) direct FDTD numerical 
simulations for the 2D PhC structure (see Fig. lU with defects separated by 2 lattice rows (point #2 in Fig. [U logio [^^/7±] — —0.481) 
for R < Re (top) and R > Re (bottom). For the Class-C coupled mode theory results, the red dashed and blue solid lines show the mode 
amplitudes |^i,2(t)|, respectively. Initially both modes are excited and the second mode is given an advantage (^1-2 =2:3). The gray 
line shows the electric field at the mode's maximum \Er(t)\ = i^max \Ei(t)e^+ + E2(t)e^- |. For the FDTD results, the field envelope at the 
center of either defect \E(rc,t)\ is shown (sampled at local extrema of light oscillations). The insets show an enlarged portion of the plots to 
show the 2Auj intermodal beating whenever both modes lase at the same time. 



cording to Eq. ^. Taking the absolute value, light oscillations 
are also neglected, so the results can be compared to the sim- 
ulations. 

In all examples of Figs. [TOl [TT] (which correspond to the 
laser operating way above threshold), the field dynamics 
shows a good qualitative and quantitative correspondence. 
Below Re where simultaneous two-mode lasing is expected, 
the in-cavity field envelope shows the characteristic 2 Acj beat 
oscillations (see the insets in Figs. fTQl - [TT1) , marking the pres- 
ence of both modes in the laser radiation. Above R^ the 
steady- state envelope is flat, indicative of single-mode lasing, 
and the beat oscillations are seen to vanish. This corresponds 
to quenching of the weaker mode in agreement with theoreti- 
cal expectations in the bistable regime. 

Some quantitative discrepancies between the model and 
simulation results can be noticed. Some of them (e.g., tempo- 
ral shifts of the spikes in Fig.fTTI) result from minor deviations 
in parameters between the real simulated structure and an ide- 
alized two-mode system considered. These deviations can be 
compensated for by fine-tuning the model 1 25] . Other discrep- 
ancies such as the the difference in the field amplitudes (both 
at spike maxima and in the steady- state) can be attributed to 
gain saturation, which may introduce correction to the form of 
the expansion ([5]) for the gain medum polarization. This is a 
limitation inherent in the present coupled-mode model. How- 
ever, Eqs. (r7l)-([T0b and ([T3l)-([T4b are clearly seen to provide 
a valid description of laser mode dynamics scenario for rela- 
tively strong pumping, unlike the near- threshold (third-order 



nonlinearity) theories which are reported to fail badly in this 
regime (see 1201 ). just like the Class- A equations ([3T]) would. 
One can overcome this limitation, e.g., following the approach 
in Refs. 1120, [21I] where a generalization of Eqs. ©-([5]) is in- 
troduced. A very good agreement with numerics is reported 
recently Oil] . However, only the time-independent (steady- 
state) theory is formulated so far. 



The knowledge that stronger pumping can restore a laser 
into the bistable regime for higher Acc; is important in the de- 
sign of a laser that can have its wavelength switched by a large 
value (such as several tens of nanometers in Refs. uji \lm) • 
A rigorous explanation of this result is yet to be given. In- 
tuitively, stronger pumping rates cause shorter lasing onset 
times compared to the cavity round-trip time, so the domi- 
nation of the stronger mode can occur before the modes have 
a chance to balance themselves through the cavity. Indeed, 
it could be noticed that the transition from simultaneous to 
bistable lasing around Re is accompanied by the disappear- 
ance of 4Acc;-pulsations in the phase of some dynamical vari- 
ables. This suggests that shorter onset due to stronger pump- 
ing allows some of the variables to become phase-locked, 
which in turn influences the whole character of the mode in- 
teraction (as was seen when the transition from Eq. dZTl) to 
Eq. ([3TI) was discussed). The detailed investigation of this ef- 
fect is a subject for further studies. 



17 




0,2 0.4 0.6 

Time [ns] 



2 
L5 

1 
0.5 



> 



^ 




-lU*. 




0.6 0.1 




0.2 0.4 0.6 

Time [ns] 



Figure 11: (Color online) Same as Fig.[Tolbut for the 2D PhC structure with defects separated by 4 rows (point #4 in Fig.[8l logio [^^/li 
-1.301). 



V. CONCLUSIONS AND OUTLOOK 

In this work, we have addressed the problem of bistabil- 
ity in a microlaser by systematically formulating the coupled- 
mode model without prior assumptions on the mode or cavity 
geometry, other than the requirement of the mode orthogo- 
nality in the cavity as well as in the gain region as described 
by Eqs. © and ([T5l) . The governing equations have been de- 
rived for all laser classes, Eqs. ©-([TOl), ([T3])-([T4|), and (|27]) 
for class-C, B, and A, respectively. The issue of classifying 
the laser dynamics in the multimode case has been revisited 
taking into account the intermode frequency separation Auj as 
a new parameter influencing the laser dynamics. The model 
has been derived for the case of two modes; however, its ex- 
tension to the case of several modes can be performed along 
the same lines. 

The simplest case of the class-A laser equations has been 
analytically investigated. It has been shown that coherent 
mode interaction processes (population pulsations) can pro- 
vide an additional mode coupling channel besides incoherent 
mode interaction (spatial hole burning). This additional cou- 
pling is what brings the laser into the bistable regime, allow- 
ing the lasing mode to be chosen on demand by the initial 
condition of the cavity. This result agrees with the early the- 
oretical predictions [6, TJ. However, microcavity modes can 
have a far better matched intensity distribution inside the gain 
region, see Eq. ([18]), compared to bulk-cavity modes, which 
are usually heavily out of match unless Auj = 0. As such, 
only a moderate amount of pulsation-induced mode coupling 
is enough to enter the bistability regime in the case of a mi- 
crolaser. This means that microlasers can be bistable in a far 
greater parameter range than bulk-cavity lasers, e.g., for much 
larger Aco (Fig. [21). We have also shown that a sizable mis- 



match in the system parameters that favors one of the modes 
can destroy any chance of bistable operation. However, a 
mismatch with respect to one parameter can be compensated 
for by a mismatch with respect to another (Fig. [5]). Again, 
due to better matched intensity distributions of microcavity 
modes the bistable regime is more tolerant to such perturba- 
tions (Fig. [7]). 



In the more general class-B or class-C laser models, we 
have shown numerically that even when Aco is too large to 
allow bistability in the near-threshold class-A case, it can be 
overcome by increasing the pumping rate (Fig. [8]). The results 
of the theory are confirmed by direct numerical FDTD sim- 
ulations and are shown to be qualitatively valid for pumping 
rates several orders of magnitudes above the lasing threshold. 
Further results on bistability in class-B/C microlasers will be 
available in a forthcoming publication. 



Bistable operation of a multimode microlaser can be use- 
ful in many respects. Since there is no need for an exter- 
nal (and potentially slow) cavity-tuning process, ultrafast all- 
optical mode switching mechanisms can be imagined. Such 
switching, occurring across c^ 20 nm on a picosecond time 
scale had indeed been demonstrated numerically in our ear- 
lier work [17]. The fast switching between stable states and 
the relatively low power of microlasers can be used in the de- 
sign of an optical memory (flip-flop) cell. We believe, that a 
compact-sized microlaser capable of multiple- wavelength op- 
eration in a wide wavelength range can find numerous appli- 
cations in integrated optics and optical communication. 
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